#ifndef GRID_HPP
#define GRID_HPP
#include "Types.hpp"
#include <cmath>
#include <algorithm>

template<typename Type, uint Cellsize >
class Grid
{
public:
    Grid();
    Grid( uint xsize, uint ysize );
    Grid(const Grid<Type,Cellsize>& s);
    ~Grid()
    {
        delete [] data_;
    }

	void reSize(uint xnew, uint ynew);
    Grid& operator=(const Grid<Type,Cellsize>& s);
    inline const uint& getxSize() const{return xsize_;}
    inline const uint& getySize() const{return ysize_;}
    inline Type& operator()( uint x, uint y, uint f );
    inline Type& operator()( uint x, uint y, uint f )const;
    void swap( Grid& grid ) {
        std::swap( ysize_, grid.ysize_ );
        std::swap( xsize_, grid.xsize_ );
        std::swap( data_, grid.data_ );
    }   
private:
    uint xsize_; // Number of nodes in x-dimension
    uint ysize_; // Number of nodes in y-dimension
    Type* data_; // Linearized, 1-dimensional representation
    // of the 2D data grid
};

// Implementation of the default constructor
template< typename Type, uint Cellsize >
Grid<Type,Cellsize>::Grid(): xsize_(0), ysize_(0), data_(new Type[1]) {}
// Implementation of the initialization constructor
template< typename Type, uint Cellsize >
Grid<Type,Cellsize>::Grid( uint xsize, uint ysize ) : xsize_(xsize), ysize_(ysize), data_( new Type[Cellsize*xsize*ysize] ) 
{
    for(uint i=0;i<Cellsize*xsize*ysize;++i)
        data_[i]=0;
}

template< typename Type, uint Cellsize >
Grid<Type,Cellsize>::Grid(const Grid<Type,Cellsize>& s):
xsize_(s.xsize_), 
ysize_(s.ysize_)
{
    delete [] data_;
    data_ = new Type[Cellsize*xsize_*ysize_];
    for(uint i=0;i<Cellsize*xsize_*ysize_;++i)
        data_[i]=0;
}

template< typename Type, uint Cellsize >
void Grid<Type,Cellsize>::reSize(uint xnew, uint ynew)
{
    delete [] data_;
    data_ = new Type[Cellsize*xnew*ynew];
    for(uint i=0;i<Cellsize*xnew*ynew;++i)
        data_[i]=0;
}

template< typename Type, uint Cellsize >
Grid<Type,Cellsize>& Grid<Type,Cellsize>::operator=(const Grid<Type,Cellsize>& s)
{
    delete data_;
    xsize_=s.xsize_;
    ysize_=s.ysize_;
    data_ = new Type[Cellsize*xsize_*ysize_];
    for(uint i=0;i<Cellsize*xsize_*ysize_;++i)
        data_[i]=0;
    return *this;
}

// Implementation of the function call operator
template< typename Type, uint Cellsize >
inline Type& Grid<Type,Cellsize>::operator()( uint x, uint y, uint f )
{
    //assert( x < xsize_ && y < ysize_ && f < Cellsize );
    return data_[y*xsize_*Cellsize+x*Cellsize+f];
}
// Implementation of the const function call operator
template< typename Type, uint Cellsize >
inline Type& Grid<Type,Cellsize>::operator()( uint x, uint y, uint f ) const
{
    //assert( x < xsize_ && y < ysize_ && f < Cellsize );
    //collision-optimized data layout
    return data_[y*xsize_*Cellsize+x*Cellsize+f];
    //propagation-optimized data layout; not testet yet
    //return data_[xsize*x_+y+xsize*ysizef];
}

// Partial template specialization for Cellsize = 0
// No class definition => compile time error
template< typename Type >
class Grid<Type,0> {};


// Partial template specialization for Cellsize = 1
template< typename Type >
class Grid<Type,1>
{
public:
    Grid();
    Grid( uint xsize, uint ysize );
    ~Grid()
    {
        delete [] data_;
    }
    Grid(const Grid<Type,1>& s);
	void reSize(uint xnew, uint ynew);
    Grid& operator=(const Grid<Type,1>& s);
    inline const uint& getxSize()const{return xsize_;}
    inline const uint& getySize()const{return ysize_;}
    inline Type& operator()( uint x, uint y );
    inline Type& operator()( uint x, uint y ) const;
private:
    uint xsize_; // Number of nodes in x-dimension
    uint ysize_; // Number of nodes in y-dimension
    Type* data_; // Linearized, 1-dimensional representation
};

// Implementation of the default constructor
template< typename Type >
Grid<Type,1>::Grid(): xsize_(0), ysize_(0), data_(new Type[1]) {}
// Implementation of the initialization constructor
template< typename Type >
Grid<Type,1>::Grid( uint xsize, uint ysize ) : xsize_(xsize), ysize_(ysize), data_( new Type[xsize*ysize] ) 
{
    for(uint i=0;i<xsize*ysize;++i)
        data_[i]=0;
}

template< typename Type>
Grid<Type,1>::Grid(const Grid<Type,1>& s):
xsize_(s.xsize_), 
ysize_(s.ysize_)
{
    delete [] data_;
    data_ = new Type[xsize_*ysize_];
    for(uint i=0;i<xsize_*ysize_;++i)
        data_[i]=0;
}


template< typename Type >
void Grid<Type,1>::reSize(uint xnew, uint ynew)
{
    delete [] data_;
    data_ = new Type[xnew*ynew];
    for(uint i=0;i<xnew*ynew;++i)
        data_[i]=0;
}


template< typename Type>
Grid<Type,1>& Grid<Type,1>::operator=(const Grid<Type,1>& s)
{
    delete data_;
    xsize_=s.xsize_;
    ysize_=s.ysize_;
    data_ = new Type[xsize_*ysize_];
    for(uint i=0;i<xsize_*ysize_;++i)
        data_[i]=0;
    return *this;
}

// Implementation of the function call operator
template< typename Type>
inline Type& Grid<Type,1>::operator()( uint x, uint y)
{
    return data_[y*xsize_+x];
}
// Implementation of the const function call operator
template< typename Type >
inline Type& Grid<Type,1>::operator()( uint x, uint y) const
{
    return data_[y*xsize_+x];
}

//global swap function
template< typename Type, size_t N >
inline void swap( Grid<Type,N>& a, Grid<Type,N>& b )
{
    a.swap( b );
}

typedef Grid<real,9> PDF_Field;
typedef Grid<real,2> V_Field;
typedef Grid<real,1> D_Field;
typedef Grid<uint,1> Flags;



#endif //GRID_HPP
